Matlab tutorial Part 2
ME 461
Professor Tilbury
Note: These tutorials have been completely revised and updated. See Control Tutorials for Matlab
for the latest version.
Representations of a linear system
There are several different ways to describe a linear system. The state-space
representation is given by the equations: .
x = Fx + Gu
y = Hx
where x is an n-vector, F is an nxn matrix, G is an n-vector, u and y are
scalars, and H is an n-row vector. (The vector x is called the state). Matlab
uses the abbreviation ss for state-space. From the previous tutorial, we
have the state-space representation: (see the tutorial
part 1 for the values of the constants Ks, Kw, b, m1, m2) F = [0 1 0 0
-(Ks+Kw)/m1 -b/m1 Ks/m1 b/m1
0 0 0 1
Ks/m2 b/m2 -Ks/m2 -b/m2];
G = [0
Kw/m1
0
0];
H = [1 0 0 0];
J = [0];Make an m-file which contains the constant values as well as these
4 matrices.
The transfer function representation is the Laplace transform of the
output divided by the Laplace transform of the input (for zero initial state):
Y(s) b_0 s^m + b_1 s^(m-1) + ... + b_m
H(s) = ------ = ------------------------------------------
U(s) s^n + a_1 s^(n-1) + ... + a_(n-1) s + a_n
The notation "^" should be read as "to the power of" and the notation "_"
should be understood as a subscript. Matlab uses the abbreviation tf for
transfer function. Since matlab can only deal with real numbers, and not symbols
(such as the Laplace variable s), a transfer function is represented in matlab
as two vectors which contain the coefficients of s in the numerator and
denominator polynomials: num = [b_0 b_1 ... b_m]
den = [a_0 a_1 a_2 ... a_n]
By convention, a_0 is always equal to 1 (if it is not, then scale the
numerator and denominator by a_0). The transfer function of the suspension
system can be found using the matlab command ss2tf
which takes the state-space representation and converts it to the transfer
function representation. Add the following line to your m-file: [num,den] = ss2tf(F,G,H,J);
and run the file. This will assign the coefficients of the numerator
polynomial to the varilable "num" and the coefficients of the denominator
polynomial to the variable "den". The command step(num,den,t)
should give exactly the same output as step(F,G,H,J,1,t)
did on the previous tutorial. Verify this now.
There are many different ways to do this verification. You can get
both plots to show up on the screen at once and compare them with your eyes.
The figure command will open up a new plot window on the screen
and make it the default for the next plot. Thus, in your m-file, you can have
the commands step(F,G,H,J,1,t);
[num,den] = ss2tf(F,G,H,J);
figure;
step(num,den,t);
and the plot which is the output of the second step
command won't overwrite the plot in the first command.
You could put both plots in the same window. The command
subplot(m,n,p)
breaks the current graph window into an mxn array of subplots and
chooses the pth one to be current (the counting is done across the top row
first, then to the second row, etc). For this example, you might try: subplot(2,1,1);
step(F,G,H,J,1,t);
title('Step response from state-space equations');
[num,den] = ss2tf(F,G,H,J);
subplot(2,1,2);
step(num,den,t);
title('Step response from transfer function');
The two responses should come out on top of each other, like: 
You could store the output in two different vectors:
yss = step(F,G,H,J,1,t);
ytf = step(num,den,t);
and then either plot them on the same graph or find the difference
between them. The difference can be found pointwise using a for loop: temp = 0;
for i = 1:length(t), temp = temp + (yss(i) - ytf(i))^2, end;
error = sqrt(temp);
Note particularly where the commas go in the for
syntax.
One of the advantages of using matlab is to not have to muck around
with the indices of a vector, but just deal with the entire vector. Try
diff = yss - ytf;
error = sqrt(sum(diff.^2));
Note that the "." before the exponent means
take the square of each element of the vector. You would get an error message
in this case if you left out the "." (try it). Just for grins, define the 2x2
matrix A = [1 2; 1 3];
and then compare the two results A^2
A.^2
Of course, even this is overkill. Matlab can take the norm of a vector
very easily: error = norm(yss - ytf)
How big is the numerical error? In theory, the two are exactly equal and
you should get 0.
Let's look at the numerator and denominator of the transfer function. To the
matlab prompt, type
den
The output should look like: den =
1.0e+06 *
0.0000 0.0010 0.0510 2.0000 2.0000
Hmmm, it looks like the first term is 0, whereas above, we said that a_0 =
1 by convention. Double-check this by displaying the first term only, den(1)
There are many different ways to have matlab display the output.
You can control this with the format command. For example, format short
is what was first shown above; it gives 5 significant digits. Because of
the magnitude of the numbers involved (>10000), matlab put a constant
factor outside of the vector. There is also format long
which displays more significant digits. When dealing with numbers of
greatly varying magnitudes, format short e
is perhaps most convenient; this gives the output: den =
1.0000e+00 1.0400e+03 5.1040e+04 2.0000e+06 2.0000e+06
If you would always like to see your numbers in this format, put the
command "format short e" in an m-file called "startup.m" in your me461
directory. Every time you start matlab from within that directory, all of the
statements in that file will automatically be executed.
Now let's
look at the numerator. The first term is zero and the second term is effectively
zero (you can start to see how numerical problems can arise when using software
like matlab).
Matlab has a variable called "eps" which is a tolerance for
determining singularity, etc. You can check the value of eps by typing eps
to matlab. This is the smallest value that matlab can handle;
anything less than this is considered to be zero. It is useful for comparing
how "close" to zero small numbers are. How close to zero is num(2)? How close
to zero was the error between the two step responses?
The numerator
polynomial could be interpreted as m = 4 and b_0 = 0, b_1 = 0, but more
traditionally, we would say that m = 2. The transfer function for this system is
therefore: 50000 s^2 + 200000 s + 2000000
H(s) = --------------------------------------------------
s^4 + 1040 s^3 + 51040 s^2 + 2000000 s + 2000000
Matlab can also find the partial fraction expansion of a transfer function
by the "residue" function. If you really need to do this you can read the help
page in matlab; pay particular attention to the warning about numerical
ill-conditioning.
Analytically, the transfer function is found from the state-space
representation using the equation: -1
H(s) = H (sI - F) G
where I is the nxn identity matrix. From this equation, it can be seen
that the denominator of H(s) is equal to the determinant of the matrix (sI -
F). Thus, if there are n states, then the denominator polynomial should have
degree n. In this case, we started with 4 states and got a denominator
polynomial with degree 4 (that is, 5 coefficients). You can check the degree
of the denominator easily by giving matlab the command length(den)
This is a good way to check for obvious mistakes.
The third
way to represent a linear system is by a list of its poles and zeros. Because
the transfer function is a rational polynomial function, its numerator and
denominator can both be factored: (s - p_1)(s - p_2) ... (s - p_m)
H(s) = k ------------------------------------
(s - z_1)(s - z_2) ... ... (s - z_n)
Note first of all that m and n here are the same as m and n from above.
Also, because of the way that we factored the two polynomials, there might be a
constant k (which will in fact be the same as b_0 which was defined above) which
appears in front of the transfer function. The transfer function is completely
defined by its list of poles and zeros together with the constant gain k. Matlab
uses the abbreviation zp to denote the pole-zero (or
zero-pole?) representation. The function tf2zp
can be used to make matlab factor the numerator and denominator
polynomials of a transfer function H(s). It returns two vectors which are the
lists of poles and zeros as well as a real number, k. Try giving matlab the
command: [z,p,k] = tf2zp(num,den)
and look at the list of poles and zeros. One of the zeros is very large
negative; this probably comes from the fact that the numerator polynomial had a
very small value for num(2) which we thought might be due to a numerical error.
Now try
[z,p,k] = ss2zp(F,G,H,J,1)
and check how many zeros there are. You might also want to see what kind
of a transfer function you can get using the function [num1,den1] = zp2tf(z,p,k}
Check the numerator in particular and compare it to the previous result.
Note that there also exist functions in matlab
tf2ss
zp2ss
which work in a hopefully analogous fashion to those described above (or
use "help" for more information).
One thing to take note of is that while the transfer function and pole-zero
representations are uniquely defined for a linear system, there are MANY
state-space representations which give the SAME transfer function! The transfer
function only encodes the input-output behavior of the system; there are many
ways to define a set of states which will achieve this behavior. We will discuss
this more towards the end of the class, in Chapter 7.
For an example, give matlab the command:
[f,g,h,j] = zp2ss(z,p,k)
with the values of z,p,k that you calculated above and look at the output.
The matrices should NOT be the same as F,G,H,J. Then try [num1,den1] = ss2tf(f,g,h,j)
and verify that this is the same transfer function that you got before.
I hope this exercise has given you an appreciation for the usefulness of matlab
as well as the usefulness of pencil-and-paper calculations and common sense to
be able to check the output of matlab!
Go to Matlab
Tutorial part 3
ME 461
Home Page
Last modified: 23 january 1995/tilbury@umich.edu
Please send any comments, clarifications, suggestions, or corrections to:
tilbury@umich.edu. Thanks! I hope you found this document helpful.